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Modeling of the subgrid-scale term of the filtered 
magnetic field transport equation 

By G. Balaracf, A. G. Kosovichev^, O. Brugieref, A. A. Wray^[ and N. N. 

Mansour^f. 

Accurate subgrid-scale turbulence models are needed to perform realistic numerical 
magnetoliydrodynamic (MHD) simulations of the subsurface flows of the Sun. To perform 
large-eddy simulations (LES) of turbulent MHD flows, three unknown terms have to be 
modeled. As a first step, this work proposes to use a priori tests to measure the accuracy 
of various models proposed to predict the SGS term appearing in the transport equation 
of the filtered magnetic field. It is proposed to evaluate the SGS model accuracy in 
term of "structural" and "functional" performance, i.e. the model capacity to locally 
approximate the unknown term and to reproduce its energetic action, respectively. From 
our tests, it appears that a mixed model based on the scale-similarity model has better 
performance. 

1. Introduction 

Significant progress towards quantitative understanding of the Sun and predictive ca- 
pabilities for solar activity and space weather requires large-scale, integrated modeling of 
the physical conditions in subsurface layers of the Sun. Realistic numerical magnetoliy- 
drodynamic (MHD) simulations of the subsurface flows and magnetic structures have be- 
come achievable because of the development of fast supercomputer systems and efficient 
parallel computer codes. These simulations are extremely important for understanding 
the complicated physics of the upper turbulent convective boundary layer of the Sun. The 
dynamics of this layer is critical for understanding of the formation of magnetic regions 
on the Sun and their activity. This layer is also a source of solar oscillations. The wave 
excitation and propagation properties change dramatically in strong field regions. Their 
modeling and investigation are very important for local helioseismology and helioseis- 
mic data analysis. Accurate Large Eddy Simulations (LES) of the subsurface dynamics 
depend on the development of specific and accurate subgrid-scale turbulence models. 
These models have to provide a realistic description of effects of small-scale unresolved 
turbulence, which are particularly important for studying wave excitation. 

In LES of MHD flows, to know the filtered velocity, Uj, and magnetic, fields, the 
filtered MHD equations, expressed in Alfven-speed units, have to be solved, 
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where D is magnetic diffusivity, v is kinematic viscosity, p is filtered total pressure. In 



these equations, 
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are subgrid-scale (SGS) tensors. In a LES, these SGS tensors cannot explicitly be deter- 
mined but are estimated via SGS models assuming relationships with resolved quantities. 
Various works have addressed the modeling of the three SGS tensors of LES of MHD fl ows 
( Yoshizawalll990nTheobald et aLl[l99l iMiiller fc Caratill2002t iMiki fc Menonll2008l ). In 
this work, we propose to investigate the performance of the models proposed for r-"- 6 from 
a priori tests. The goal is to distinguish the modeling error due to the models of rff from 
the other modeling errors. The modeling performances are evaluated as "structural" and 
"functional" performances. 
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2. Modeling of rf^ and performance measurement 

2.1. Available SGS models 

The modeling of the unknown SGS term, rfj , appearing in the filtered transport equation 
of the magnetic field has been addressed in various ways in the past. First, an approach 
based on the definition of an eddy magnetic diffusivity, D t , has been proposed, leading 
to the general gradient-diffusion model expression 



= -2 A J, 



(2.1) 



with Jij = i ( — J^-J , the filtered magnetic rotation tensor. Various definitions of D t 

are a vailable in the literature. Va rious extens i ons of the Smagorinsky model ( Smagorinskvl 
1963[) have been proposed. First, Yoshizawa ( 1987 ) defined the eddy magnetic diffusivity 

as 

_ \ 1/2 

Cxjf , (2.2) 
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D t = C\A 2 y-CvSl 

\ (fe" + if 1 )' tne resorve d velocity rate-of-strain tensor and j = V x 6, 
the current density. Assum ing a local equilibrium between production and dissipation, 
Hamba fc Tsuchival (|2010T) determine the constant values as C\ — jC v and C v = 0.046. 
Theobald et all (|1994) propose to define D t only with the current density norm, 



with Sij — 



D t = dA 2 \j\. 



(2.3) 



Miiller fc Garni! (2002) use this model by computing d ynamically the model coefficien t 
as usually done fo r LES of hydrodynamic turbulence! Lilly! 1 1992 : Germano et al 1 ll99lh . 
In the same paper, iMuller fc Caratl (|2002h define a new eddy magnetic diffusivity based 
on the cross helicity dissipation, 



=c 2 A 2 ij.<3r, 



(2.4) 



with to the vorticity vector. The model coefficient is also computed dynamically. 

Another approach to model is based on the filtering operation itself. For example, a 
Taylor series expansion of a filtered product, fg (where / and g are quantities des c ribin, 
magnetic or flow fields), can be given for a Gaussian filter as (see iBalarac et al\ ( 200. 
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Keeping only the first term, a gradient model for r z " h can be written as 
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Another model can be constructed from the scale-similarity hypothesis proposed by 
( 19801) . The main idea is to assume that the statistical structure of the 
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tensors constructed on the basis of the subgrid scales is similar to that of their equiv- 
alents evaluated on the basis of the smallest resolved scales. From this assumption, the 



scale-similarity model for rff can be written as 
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However, these types of models are known to not provide enough energy transfer between 
grid scale (GS) and subgrid-scale (SGS), leading t o unstable simulation s when they is 
used to close the filtered Navier-Stokes equations (jVreman et "all ll997l). To avoid this 
unstable behavior, a mixed model can be proposed. As proposed by Clark et all ( 19791 ) 
for the Navier-Stokes equations, the mixed model consists to add an eddy magnetic 
diffusivity approach to the gradient or the scale-similarity model. Using the model (|2.3j) . 
the mixed-gradient model is thus written as 
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whereas the mixed-scale-similarity model is written as 
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with C3 and C4 dynamically computed. 

Fin ally, a last approach i s inspired from works in RANS context ( Yoshizawa 1990f ). 
Thus. lMiki fc Menonl (|2008h proposed to model r\f as 



(2.10) 



where ekij is the Levi-Civita symbol and Ek is the turbulence electromotive force defined 
as 

E k = ab k - (3j k + juj k . (2-11) 

In this equation, a, (3 and 7 are spatially dependent quantities defined with the kinetic 
and magnetic subgrid-scale energy, k sgs = \i^Ui~UiUi) an d fc sgs,fc = ^ (b j bj — b jbj) , 
and computing model coefficient locally and dynamically (see Miki fc Menon ( 20081) for 
details). Note that this model needs to solve two additional transport equation for k sgs 
and k S9S,b with additional unknown terms to model. 



Eight SGS models for rS presented above have been summarized in table 12.11 Some 
of these models are very easy to apply (especially without dynamic coefficient computa- 
tion) but some others need a non-negligible implement effort and have a non-negligible 
computational cost. The goal of this work is to be able to determine which model appears 
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Name 

Yoshizawa 
Theolbald 
Cross-helicity 

Gradient 
Scale-similarity 
Mixed-Gradient 
Mixed-Scale-similarity 
Miki 



Equations 



1 23)) and (2T2 1 



( 2.1 1 and (2,31 



( 2TT1) and ((274 1 
11231) 
(I2T1) 

im 

(pAO) and (l2~TTj) 



Symbol 



□ 



A 

V 

A (dotted line) 
V (dotted line) 
> 



Remark 

2 constant coefficients 
1 dynamic coefficient 
1 dynamic coefficient 
1 constant coefficient 

Explicit filtering 
1 dynamic coefficient 
1 dynamic coefficient 
2 additional transport equations 



Table 1. Models for the unknown SGS term rf. 



as the best way to model T 4 " b . Our analysis is mainly based on a priori tests using di- 
rect numerical simulation (DNS) data. Our a priori tests are performed to measure two 
types of per formance. Usi ng the distinction between structural and functional models 
proposed bv lSagautl (|2005I l we propose to measure the SGS model performance in terms 



of structural and functional performance. 

2.2. Structural performance and optimal estimator 

The structural modeling strategy consists in making the best approximation of the un- 
known SGS term by constructing it from the knowledge of the structure of small-scales. 
For example, the gradient model, Eq. (12.6[) . is a structural model based on a Taylor series 
expansion of the filtering operation. This type of model is known to give a good approx- 
imation of the unknown term with a high correlation between the unknown term and 
the model in a priori test. From this definition, we define the structural performance of 
a model by its capacity to locally approximate the SGS term to be modeled. 

To evaluate the structural performanc e, the optimal estimation theory is used. In the 
framework of optimal estimation t heory (|Deutsch[ [l965). the models are compared using 
the notion of an optimal estimator (jMoreau et al\\2Q Q&\ . Based on this idea, if a quantity 
t is modeled with a finite set of variables <fi, an exact model cannot be guaranteed. If the 
exact solution r is known, for example from DNS, the optimal estimator of r in terms 
of the set of variables <p is given by the expectation of the quantity r conditioned on the 
variables in the set, i.e. (r|</>), where the angular brackets indicate statistical averaging 
over a suitable ensemble. A quadratic error can consequently be defined as the average 
of the square of the difference at each point between the conditional mean value given 
by the value of <p at this point and the exact value of the quantity, 

emi„ = ((r-(r|0)) 2 ), (2.12) 

where e m i n is the error. It should be noted that any model formulated using the variable 
set <j) will introduce an error that is larger than or equal to this minimum error, with 
the best model formulation producing this minimum error. Consequently, this quadratic 
error e min is referred to as the irreducible error. Only a change in the variable set may 
reduce the magnitude of this error. In contrast, the total quadratic error is given as 

etot = ((r-/(0)) 2 ), (2.13) 

with f((f>) the proposed model for r. The method allows to compare different LES models 
by comparing their total errors. The method allows also to evaluate the improvement 
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Figure 1. Kinetic, E(k), and magnetic, Es(k), energy spectra: E(k) ; Esik) • 

The thin vertical lines show the location of some filters used in this work: A/ Ax —2, 4, 8 and 
16. 

possibility of a given model without changing of the set of parameters. Indeed, if the 
total error of a given model is much higher than its irreducible part, improvement can 
be expected (by modification of the coefficient computation, for example). This method 
finally allows to find the best set of quantities to model a sub-filter term by comparing 
their irreducible error. 

2.3. Functional performance 

The functional modeling strategy consists to consider the action of the subgrid terms on 
the transported quantity (here, the magnetic field) and not the unknown term itself. It 
can consist in introducing a dissipative term, for example, that as a similar effect but 
not necessarily the same structure. Since an adequate mechanism to dissipate magnetic 
energy from resolved to subgrid scales is essential, we define the functional performance 
as the model capacity to lead to the good energy dissipation. The transport equation of 
the GS magnetic energy, k b = \bjbi, is 

dk b dk b 1 _ d 2 k b _ dbi dh d - _ , d ,-rs d .- ub . ub - 

In this equation, the GS/SGS energy exchanges is due to the SGS dissipation, r^Jij. 
The functional performance will be then evaluated as the model capacity to well predict 
this quantity in a statistical sense. 

3. Results 

3.1. Numerical method 

As already explained, to clearly identify the performances of the various modeling strat- 
egy of T™ b without taking into account the modeling error of the other SGS unknown 
terms appearing in the filtered MHD equations (11.1[) - (|1.3[) . we performed a priori tests. A 
priori tests are conducted using direct numerical simulation (DNS) data from a forced ho- 
mogeneous isotropic turbulence computation. A pseudo-spectral code with second-order 
explicit Runge-Kutta time-advancement is used. The viscous terms are treated exactly. 
The simulation domain is discretized using 256 3 grid points on a domain of length 2ir. 
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Figure 2. Structural performance of models. Evolution of the irreducible and total errors, 
Eq. (|2.12|) and (|2.13|) . with the filter size for box (left) and spectral cutoff (right) filter. Open 
and solid symbols are for total and irreducible errors, respectively (see Table 12.11 for symbols 
correspondence) . 



A classic 3/2 rule is used for dealiasing the non-linear conv ection term, and statistical 
stationarity is achieved using a forcing term ( Alvelius 19991 ). The transport equation of 



the magnetic field is advanced simultaneously using an identical numerical scheme. First, 
a hydrodynamic (no magnetic field) case is performed; when the statistically stationary 
state is obtained, the Reynolds number based on the Taylor microscale is around 100. 
The magnetic field is then initialized at small scales with a small amplitude. The mag- 
netic Prandtl number is set to 0.5. Without external forcing, the magnetic energy grows 
leading quickly to a new statistical state. The a priori tests are performed when the flow 
is statistically stationary. 

In the a priori tests, explicit filter is used to replicate the behavior of the filter implicitly 
associated with the discretization in real LES. Two kinds of filter are used, a spectral 
cutoff filter to mimic spectral LE S and a box filter to mimic LES using centered finite 
differences (jRogallo fc Moinl 1984} . Several different filter sizes have been used such as 



2 < A/Ax < 16, where A is the filter width and Ax is the grid spacing used in the DNS. 
The location in wavenumber space of the filters used are displayed in Fig. [I] which shows 
the kinetic and magnetic energy spectra. 

3.2. Models performances 

As already explained, to first evaluate the structural performance of the SGS models, 
the total error, Eq. (I2.13[) for each model is considered and compared with its irreducible 
error, Eq. (12.121) . Figure [2] shows the evolution of the total and irreducible error with 
the filter size, for the different SGS models and for both box and spectral cutoff filter. 
On this figure, the error is normalized by the statistical variance of the SGS term. First 
conclusions can be addressed. As expected, SGS models based on structural approach 
(models based on scale-similarity assumption or on Taylor series expansion) lead to the 
smallest errors. Note that for both filters the models based on the scale-similarity as- 
sumption have the smallest errors. In particular for the spectral cutoff filter the errors 
of the models based on a Taylor series expansion stay high. This is becau se the spectra l 
cutoff filter leads to a divergent Taylor series, because of its non-localness (|Sagautll2005l ). 
All the other models have irreducible errors higher than the total error model of the 
models based on structural approach. This shows that a structural improvement of these 
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Figure 3. Functional performance of models. Evolution of the mean GS/SGS magnetic energy 
transfer with the filter size for box (left) and spectral cutoff (right) filter (see Table 12.11 for 
symbols correspondence). The thick line is the mean GS/SGS magnetic energy transfer from 
the filtered DNS. 



models can not be expected without adding new quantities in its set of parameters. This 
allows to illustrate the improvement due to the mixed approach in comparison with the 
models based only on an eddy diffusivity assumption. 

As second step, the functional performance is now studied from the evolution of the 
mean GS/SGS magnetic energy transfer, {rffJij), with the filter size. Figure [3] shows the 
results for box and spectral cutoff filters. First, it is shown that the "pure" structural mod- 
els, i.e. the gradient and scale-similarity models, give not enough transfer in comparison 
with the DNS results. This is a well-known problem for models based only on structural 
approach for hydrodynamic LES. Indeed this model as known to not give enough dis- 
sipation, leading to unstable simulations. Conversely, Theobald's eddy-diffusivity based 
model predicts enough GS/SGS transfer and even an over-estimation for the spectral 
cutoff filter in comparison with the DNS results. Thus, this allows the mixed models 
to predict enough dissipation. Note that the other eddy-diffusivity based models do not 
lead to enough GS/SGS transfer to be a good candidate to b uild a mixed model. In 
particular, the cross- helicity model leads to no GS/SGS transfer. iMuller fc Carati ( 2002 1 
had already shown this property explaining that in LES performed with this SGS model, 
the lack of GS/SGS magnetic energy transfer is compensated by the transfer between 
kinetic and magnetic energy caused by the Lorentz force. 

From this analysis, the best performing model appears to be the mixed-scale-similarity 
model. First, for the structural performance, the analysis based on the total and irre- 
ducible errors shows that we can not expect an improvement of the other models to 
have the same performance of this model. Indeed the total error of this model is always 
smaller than the irreducible error, i.e. the smallest possible error, of the other models. 
For the functional performance, the eddy-diffusivity part of this model allows to predict 
enough dissipation to allow the simulation stability. It can be however noted that a pos- 
sible improvement of this model would be to correct the over-estimation of the GS/SGS 
transfer. 



3.3. Dynamic procedure at the tensor divergence level 

In the results above, the mixed-scale-similarity model, Eq. (|2.9|) . used a classic dynamic 
procedure to compute C4. The classic dynamic procedure uses a second (test) filter, de- 
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Figure 4. Comparison of model performance for the mixed-scale-similarity model with the 
classical (triangle) and with the divergence based (cross) dynamic procedure. Left: evolution of 
the irreducible (solid line) and total errors with the filter size. Right: Evolution of the mean 
GS/SGS magnetic energy transfer with the filter size and comparison with the mean GS/SGS 
magnetic energy transfer from the filtered DNS (solid line). 



noted ?, of size A = 2 A. This procedure is based on the Germano identity (iGermano et al 



19911 ). 



r ub 



where T"f = 



(uibj - UibjJ. 



rpub 

ij 

Thus, 



T UC 



= h 



can 



be computed directly from the resolved field. Assuming that, 7$ , which is the subgrid 
tensor corresponding to the second filtering level, is also modeled with the mixed-scale- 
similarity model and with the same value of GV An equation for C4 can then be written 



from a least squares averaging procedure (|Liilvlll992D . In fact, from the filtered transport 
equation of the magnetic field, Eq. fj 1 . 2[) . it can be noted that only the vector given by 
the divergenc e of the tensor, dr ]f /dxj has to be known and not the tensor, r z " & , itself. 
In this sense, Clark et al. ( 19791 ) explained the efficiency of the Smagorinsky model for 
the Navier-Stokes equations. It is shown that the correlation between the Smagorinsky 
model and the SGS term is weak at the tensor level but higher at the vector level. Thus 
to improve the prediction of the GS/SGS magnetic energy transfer with the mixed-scale- 
similarity model, a dynamic computation of the C4 coefficient is tested for dr^ /dxj 
instead of r-f. The starting point is thus to use the divergence of the Germano identity 



U ij 



dr ub 
13 



dxj dxj dxj 

and the same steps of the classical dynamic procedure are then used. Figure 2] shows 
the structural and functional performance of the mixed-scale-similarity model using a 
divergence based dynamic procedure. The results are compared with the classic mixed- 
scale-similarity model. It is first important to note that the structural performance has 
not deteriorated. The total error of the model stays close to its irreducible error and the 
error is much smaller than the error of the other models (Fig. [2]). Moreover, the functional 
performance is improved. The new dynamic procedure allows a better prediction of the 
GS/SGS magnetic energy transfer. The over-prediction observed with the classic dynamic 
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Figure 5. Time evolution of the filtered kinetic (left) and magnetic (right) energy. DNS (solid 
line), LES with mixed scale similarity model using a classical dynamic procedure (dashed line) 
and LES with mixed scale similarity model using a divergence based dynamic procedure (dotted 
line) 



procedure disappears and the transfer predicted by the model is close to the transfer 
computed from the filtered DNS. 

Direct and large-eddy simulations have been performed for decaying MHD turbulence. 
LES's use mixed-scale-similarity models both with classical and divergence-based dy- 
namic procedures. For these LES, the computational domain is discretized using 64 3 grid 
points. The filtered Navier-Stokes equation is closed by a dynamic Smagorinsky model to 
evaluate — rA- • The influe nce of the SGS Lorentz force is just taken into account in the 
dynamic coefficient (e.g. see Miiller fc Caratil ( 2002t )). Figure [5] shows the decaying of the 
filtered kinetic (left) and filtered magnetic (right) energy for LES and filtered DNS. The 
results seem to show a better agreement with the divergence-based dynamic model, but 
it is difficult to be sure due to the coupling of the modeling error and an over-dissipation 
of the kinetic energy in the first stage of the simulations. In further work, the modeling 
improvement of r" and r 4 & - will be addressed again using a priori tests as starting point. 



4. Conclusions 

In this work, the modeling of the subgrid-scale (SGS) term appearing in the transport 
equation of the filtered magnetic field is addressed. From a priori tests, the performances 
of several SGS models have been evaluated. The measure of model performance is de- 
fined as structural and functional performances. The structural performance is defined as 
the model capacity to locally reproduce the unknown SGS term whereas the functional 
performance is defined as the model capacity to reproduce the energetic action of the 
unknown term. The structural performance is thus evaluated by using the optimal esti- 
mation theory. This allows to compare the models but also to evaluate the improvement 
possibility of a given model. The functional performance is evaluated by the comparison 
of the GS/SGS magnetic energy transfer given by the model with the expected GS/SGS 
magnetic energy transfer from the DNS data. In this work, the mixed model based on the 
scale-similarity model with a divergence- based dynamic procedure has the best perfor- 
mance. This work could be the starting point of a methodology to improve SGS modeling 
in various configurations. In further work, the modeling of other SGS terms appearing 
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in the equations of MHD flows will be addressed. 

The authors have benefited from fruitful discussions with CTR Summer Program par- 
ticipants. Computing resources were provided by IDRIS-CNRS (http://www.idris.fr/). 
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